function floating_point_demo()
% FLOATING_POINT_DEMO 浮点运算演示
% 
% 演示浮点数的表示、精度限制和舍入误差

fprintf('=== 浮点运算演示 ===\n\n');

%% 1. 机器精度
fprintf('1. 机器精度\n');
eps_machine = eps;
fprintf('机器精度 eps = %e\n', eps_machine);
fprintf('这意味着 1 + eps > 1，但 1 + eps/2 = 1\n');

% 验证
test1 = 1 + eps_machine;
test2 = 1 + eps_machine/2;
fprintf('1 + eps = %.16f\n', test1);
fprintf('1 + eps/2 = %.16f\n', test2);
fprintf('1 + eps > 1: %s\n', mat2str(test1 > 1));
fprintf('1 + eps/2 == 1: %s\n', mat2str(test2 == 1));

%% 2. 浮点数的表示范围
fprintf('\n2. 浮点数表示范围\n');
fprintf('最小正数 realmin = %e\n', realmin);
fprintf('最大正数 realmax = %e\n', realmax);
fprintf('正无穷 inf = %f\n', inf);
fprintf('非数值 NaN = %f\n', NaN);

%% 3. 舍入误差示例
fprintf('\n3. 舍入误差示例\n');

% 示例1: 大数加小数
fprintf('示例1: 大数加小数\n');
large = 1e16;
small = 1;
result1 = large + small - large;
result2 = small;
fprintf('(1e16 + 1) - 1e16 = %g\n', result1);
fprintf('理论值应该是: %g\n', result2);
fprintf('相对误差: %e\n', abs(result1 - result2) / result2);

% 示例2: 相近数相减
fprintf('\n示例2: 相近数相减\n');
x = 1.0000001;
y = 1.0000000;
diff_computed = x - y;
diff_exact = 1e-7;
fprintf('x = %.10f\n', x);
fprintf('y = %.10f\n', y);
fprintf('计算的差值: %e\n', diff_computed);
fprintf('精确差值: %e\n', diff_exact);
fprintf('相对误差: %e\n', abs(diff_computed - diff_exact) / diff_exact);

%% 4. 求和顺序的影响
fprintf('\n4. 求和顺序的影响\n');
n = 1000000;
terms = ones(n, 1) / n;  % 每项都是 1/n

% 正向求和
sum_forward = 0;
for i = 1:n
    sum_forward = sum_forward + terms(i);
end

% 反向求和
sum_backward = 0;
for i = n:-1:1
    sum_backward = sum_backward + terms(i);
end

% MATLAB内置求和
sum_builtin = sum(terms);

fprintf('正向求和结果: %.15f\n', sum_forward);
fprintf('反向求和结果: %.15f\n', sum_backward);
fprintf('内置求和结果: %.15f\n', sum_builtin);
fprintf('理论值: 1.000000000000000\n');

%% 5. 数值稳定性示例
fprintf('\n5. 数值稳定性示例\n');
fprintf('计算二次方程 x^2 - 2px + 1 = 0 的根\n');
p = 1e8;
fprintf('p = %e\n', p);

% 标准公式（数值不稳定）
discriminant = sqrt(4*p^2 - 4);
x1_unstable = (2*p + discriminant) / 2;
x2_unstable = (2*p - discriminant) / 2;

% 数值稳定的方法
x1_stable = p + sqrt(p^2 - 1);
x2_stable = 1 / x1_stable;  % 利用韦达定理

fprintf('不稳定方法:\n');
fprintf('x1 = %e\n', x1_unstable);
fprintf('x2 = %e\n', x2_unstable);
fprintf('x1 * x2 = %e (应该等于1)\n', x1_unstable * x2_unstable);

fprintf('稳定方法:\n');
fprintf('x1 = %e\n', x1_stable);
fprintf('x2 = %e\n', x2_stable);
fprintf('x1 * x2 = %e (应该等于1)\n', x1_stable * x2_stable);

end